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HIGHLIGHTS 


►  Novel  approach  to  model  solid-oxide  fuel  cell  (SOFC)  stacks  in  two  and  three  dimensions. 

►  Hierarchical  model  developed. 

►  Results  presented  for  a  stack  operated  on  pre-reformed  hydrocarbon  fuel. 

►  Transient  response  to  load  changes  studied. 
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This  paper  presents  a  novel  approach  to  model  the  transient  behavior  of  solid-oxide  fuel  cell  (SOFC) 
stacks  in  two  and  three  dimensions.  A  hierarchical  model  is  developed  by  decoupling  the  temperature  of 
the  solid  phase  from  the  fluid  phase.  The  solution  of  the  temperature  field  is  considered  as  an  elliptic 
problem,  while  each  channel  within  the  stack  is  modeled  as  a  marching  problem.  This  paper  presents  the 
numerical  model  and  cluster  algorithm  for  coupling  between  the  solid  phase  and  fluid  phase.  For 
demonstration  purposes,  results  are  presented  for  a  stack  operated  on  pre-reformed  hydrocarbon  fuel. 
Transient  response  to  load  changes  is  studied  by  introducing  step  changes  in  cell  potential  and  current. 
Furthermore,  the  effect  of  boundary  conditions  and  stack  materials  on  response  time  and  internal 
temperature  distribution  is  investigated. 

©  2012  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Solid-oxide  fuel  cells  (SOFCs)  have  been  proposed  as  a  power 
source  for  distributed  or  stationary  power  plants  and  mobile 
applications.  This  wide  usability  is  due  to  their  modularity,  high 
electrical  efficiency,  low  environmental  pollution,  fuel  flexibility, 
and  tolerance  to  fuel  pollutants.  A  typical  unit  cell  in  a  planar  SOFC 
stack  is  composed  of  a  positive  electrode-electrolyte-negative 
electrode  (PEN)  assembly,  interconnected  plates  on  both  sides, 
and  gas  seals.  In  practical  applications  of  SOFCs,  multiple  cells  are 
assembled  to  form  a  stack  and  a  serial  connection  in  an  electric  loop 
to  generate  high  voltage  and  power. 

An  attractive  feature  is  the  possibility  of  SOFC  being  combined 
with  other  power  generation  systems  (e.g.,  gas  turbines)  to  achieve 
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high  overall  electrical  power  generation  efficiency,  due  to  their  high 
working  temperatures  [1-9].  One  such  application  domain  that 
draws  considerable  interest  from  the  automobile  industry  is  the  use 
of  SOFCs  for  auxiliary  power  units  (APUs)  because  of  their  high 
efficiency  and  low  emissions.  For  SOFC-based  APUs,  it  is  not  trivial 
to  operate  the  cells  on  pure  hydrogen.  Instead,  operating  the  cells 
on  reformate  fuels  from  gasoline  or  diesel  will  be  more  attractive 
because  these  fuels  are  available  on-board  and  can  be  reformed 
before  usage  as  fuel.  Furthermore,  the  APUs  are  subjected  to 
frequent  load  changes  which  cannot  be  pre-determined.  The 
control  strategy  and  the  power  electronics  associated  with  these 
applications  may  be  quite  complex  [10-20].  Understanding  the 
response  time  of  a  system  is  particularly  important  because  of  the 
frequently  fluctuating  nature  of  demand  for  electrical  power. 

In  this  paper,  we  present  the  development  of  a  new  stack  model 
and  its  solution  algorithm.  The  model  is  demonstrated  with 
a  pragmatic  fuel  composition  for  APUs.  A  number  of  issues,  such  as 
(i)  required  power  rating  of  the  APU,  (ii)  start-up  time,  (iii)  response 
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to  load  changes,  (iv)  control  strategies,  (v)  fuel  choice,  (vi)  fuel 
utilization,  (vii)  stack  degradation,  and  (vii)  stack  volume,  need  to 
be  considered  when  addressing  the  application  of  SOFC  as  APUs. 
However,  this  paper  -  as  a  first  step  -  focuses  only  on  the  devel¬ 
opment  of  the  numerical  algorithm,  stack  simulation,  and  transient 
response  to  load  changes. 

There  is  a  great  deal  of  literature  addressing  the  dynamics  of 
load  changes  during  SOFC  operation.  Most  of  these  modeling 
studies  are  either  done  for  single  channels  [21-24],  single  cells 
[25-29]  or  for  stacks  containing  small  number  of  cells  [30-35]. 
Solving  the  electrochemical  model  equations  requires  fixing  either 
the  operating  cell  potential  or  the  current.  Li  et  al.  have  studied  the 
dynamic  response  of  the  cell  to  changes  in  fuel  flow  rates  [32].  Even 
though  Achenbach  and  Li  et  al.  simulated  cross-flow  cells,  the  time 
constant  reported  by  Li  et  al.  is  in  the  order  of  hours  and  that 
reported  by  Achenbach  is  in  the  order  of  a  few  minutes.  Achenbach 
[36,37]  examined  the  transient  cell  voltage  performance  due  to 
temperature  changes  and  perturbations  in  current  density. 
Furthermore,  Li  et  al.  also  reported  that  the  time  constants  for 
different  parts  of  the  cell  to  reach  steady  state  are  different.  Though 
both  of  these  studies  are  done  for  internal  reforming,  the  fuel 
compositions  used  differ  slightly  from  each  other.  Nevertheless,  the 
time  constant  reported  by  Li  et  al.  differs  largely  from  that  reported 
by  Achenbach.  Therefore,  it  must  be  assumed  that  the  time 
constant  for  the  cell  processes  depends  largely  on  the  model 
parameters  used.  Unfortunately,  it  is  not  possible  to  make  a  direct 
comparison  as  the  literature  does  not  provide  a  comprehensive  list 
of  material  properties  and  electrochemical  model  parameters. 

A  Laplace  domain  model  has  been  considered  by  Padulles  et  al. 
[38]  for  the  dynamics  of  a  SOFC  plant.  However,  the  model  does  not 
consider  changes  in  the  transport  fields  along  the  gas  flow  channels 
and  inside  the  electrodes.  The  model  contemplated  by  Jurado  [39] 
is  similar.  The  author  studied  fuel  cell  response  characteristics,  in 
terms  of  real  and  reactive  power,  to  frequency  step  and  voltage  step 
transients,  respectively.  The  time  constant  of  the  cell  voltage 
approximated  to  0.03-0.05  s,  although  the  exact  operating  condi¬ 
tions  of  the  cell  could  not  be  determined.  Lu  et  al.  [40]  have 
observed  that  there  can  be  differences  in  time  constants  between 
different  transport  fields.  However,  the  time  constants  of  a  partic¬ 
ular  transport  field  may  also  vary  with  different  operating  condi¬ 
tions.  Such  a  study  has  been  presented  by  Murshed  et  al.  [41  ],  in 
which  load  changes  in  current  at  different  temperatures  are 
introduced  and  the  differences  in  voltage  transients  were  observed. 
For  a  step  change  in  current,  the  time  constant  of  the  voltage 
transient  in  their  study  was  about  10-20  s,  depending  on  the 
operating  temperature.  It  may  be  mentioned  that  the  time  constant 
of  an  output  for  a  stack  is  typically  higher  than  that  for  a  single  cell. 
In  case  of  a  step  change  in  the  flow  rate  of  H2  or  02,  the  time 
constant  of  voltage  is  about  5-10  s.  For  a  step  up  in  current  density, 
the  time  constant  is  about  110-120  s,  while  in  the  case  of  a  step 
down,  it  is  about  200  s.  Iora  et  al.  [42]  have  imposed  a  ramp  load 
change  of  ±2000  A  m-2  over  a  period  of  60  s  to  study  the  transients 
in  voltage. 

The  latter  three  dynamic  stack  models  (Lu  et  al.  [40];  Murshed 
et  al.  [41];  Iora  et  al.  [42])  are  non-isothermal  and  yet  the  time 
constants  in  these  studies  are  in  seconds.  A  much  larger  time 
constant  is  observed  in  the  study  of  Hall  and  Colclaser  [43].  In  this 
non-isothermal  model,  the  time  constant  of  the  cell  terminal  voltage 
is  about  800  s  for  a  step  change  in  current  density  of  about  43%.  A  still 
higher  value  of  time  constant  is  reported  by  Lin  and  Hong  [44].  In 
response  to  a  step  change  in  current  from  0.5  to  0.3  A  cm-2,  the  time 
constant  of  the  voltage  transient  for  a  single  cell  was  found  to  be 
about  30  min.  Kazempoor  et  al.  [53]  have  investigated  the  dynamic 
characteristics  of  SOFC  systems  in  order  to  develop  relevant  control 
strategies.  They  applied  two  control  strategies,  i.e.,  cell  constant  fuel 


flow  rate  and  constant  fuel  utilization,  to  conclude  that  the  latter 
strategy  had  smaller  relaxation  times  and  temperature  gradients, 
and  was  better  suited  to  control  an  SOFC  operation.  Moreover,  they 
also  considered  an  additional  strategy  involving  introduction  of  step 
changes  in  inlet  flow  temperature,  which  resulted  in  relaxation 
times  that  were  two  times  higher. 

In  none  of  the  dynamic  models  discussed  previously,  detailed 
dynamics  of  the  diffusion  processes  inside  the  porous  electrodes 
were  considered.  However,  the  transients  are  greatly  affected  by 
the  mass  transfer  resistances  inside  the  electrodes,  as  can  be  seen  in 
the  study  of  Qi  et  al.  [26,27].  In  their  earlier  work  [26],  the  authors 
had  focused  on  cell  impedance  and  diffusion  resistances  inside  the 
electrodes.  The  diffusion  considered  is  one-dimensional  and  Tay¬ 
lor’s  expansion  is  applied  to  derive  the  flux  in  the  form  of  a  transfer 
function.  Transients  in  voltage  and  current  are  observed  in 
response  to  a  step  change  in  the  resistive  load  of  the  system.  The 
authors  concluded  that  the  slow  response  following  the  fast 
response  is  due  to  the  diffusion  limitation  inside  the  porous  layer. 
Their  subsequent  model  [27]  is  non-isothermal,  but  the  diffusion 
considered  is  still  one-dimensional  and  the  approximations  in 
deriving  the  transfer  function  model  remain  the  same  as  before. 
None  of  the  studies  presented  above  are  validated  with  experi¬ 
mental  data. 

In  the  present  analysis,  we  report  the  modeling  of  a  30  cell  stack 
made  up  of  co-flow  configuration.  The  fuel  composition  employed 
is  31.81%  H2, 36.56%  N2, 13.19%  CO,  13.19%  CH4, 2.23%  C02,  and  3.02% 
H20.  This  fuel  composition  essentially  results  from  mixing  3%  H20 
to  the  fuel  from  the  partial  oxidation  of  methane  with  air  [45]. 


2.  Modeling  approach 

2.1.  Single  unit  cell 


A  quasi-two-dimensional  model  reported  previously  is  used  for 
the  simulation  of  anode  channels,  i.e.,  the  planar  single  channel 
model  is  not  a  full  2D  model.  The  species  composition  is  resolved  in 
ID  in  the  gas  channels  (along  the  direction  of  axial  flow)  as  well  as 
in  the  porous  media  (transverse  to  the  direction  of  axial  flow).  For 
each  axial  position  in  the  flow  channel,  the  porous  media  is 
resolved  across  its  thickness,  i.e.,  even  though  the  governing 
equations  have  only  one  independent  space  coordinate,  it  gener¬ 
ates  an  overall  2D  effect.  In  this  model,  the  flow  through  the 
channels  is  assumed  to  be  plug  flow,  and  multi-component  species 
transport  through  the  porous  media  is  represented  by  Dusty  Gas 
Model  (DGM).  Detailed  description  of  the  single  channel  model  is 
also  reported  elsewhere  [46]. 


2.1.1.  Fluid  and  energy  transport 

Flow  through  fuel  and  air  channels  is  assumed  to  be  one¬ 
dimensional  and  laminar  in  nature.  The  plug  flow  equation  for 
species  continuity  in  the  channels  is  given  by 


d(pfYk) 

at 


>(pfvYk) 

6z 


+  fjkWk,  k  *1 . Kg. 


(1) 


The  velocity  is  calculated  from  the  momentum  equation 


(2) 


where  Pe  is  the  perimeter  associated  with  the  membrane  electrode 
assembly,  pf  is  the  fluid  density,  Y/<  is  the  species  mass  fraction  of 
species  k,  u  is  the  velocity,  z  is  the  axial  position,  t  is  the  time,  V\4  is 
the  species  molecular  weight,  and  Ac  is  the  cross-sectional  area  of 
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the  channel.  Assuming  constant  pressure  in  the  channels,  the 
density  is  calculated  from  the  ideal  gas  equation 


pM  =  pfRT.  (3) 

In  Eqs.  (1)  and  (2),  is  the  flux  at  the  electrode  channel  inter¬ 
face,  which  is  calculated  using  the  dusty  gas  model  (DGM)  as 
described  below.  The  single  channel  model  evaluates  the  species 
properties  based  on  the  stack  temperature  at  the  corresponding 
spatial  positions.  The  temperature  Tf  in  the  flow  channels  is 
determined  from  the  heat  balance  equation 


a 


(PfCpfTf) 

at 


a(upfcpfrf 

9z 


4h( 

Dh  {  stack 


Tf). 


(4) 


Here,  CPf  is  the  specific  heat  capacity  of  the  fluid.  Tf  and  rstack  are 
the  temperatures  of  the  fluid  stream  and  solid  phase,  respectively. 
The  first  term  on  the  right-hand  side  of  Eq.  (4)  represents  the 
transport  of  energy  due  to  bulk  fluid  flow,  and  the  second  term 
represents  the  heat  transfer  between  the  channels  and  the  solid 
stack  structure.  The  heat  transfer  coefficient  h ,  which  is  a  local 
property  and  changes  along  the  length  of  the  channel,  is  evaluated 
from  the  Nusselt  number 


(5) 


where  Dh  is  the  hydraulic  diameter  and  A  is  the  thermal  conduc¬ 
tivity  of  the  fluid.  The  Nusselt  number  is  expressed  empirically  as 
[47]: 


Nu  =  3.095  +  8.933 


1000\  -°  5386 


Gz 
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exp  - 
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where  the  Graetz  number,  Gz,  is  given  in  terms  of  the  Reynolds 
number,  Re,  and  Prandtl  number,  Pr,  as 


Gz  =  —Re  Pr. 
z 


(7) 


where  z  is  the  axial  position  along  the  channel. 


2.1.2.  Porous  media  transport 

Due  to  geometrical  considerations,  species  transport  through 
the  porous  media  is  assumed  to  be  one-dimensional  along  the 
thickness  of  the  porous  structure,  i.e.,  transverse  to  the  direction  of 
flow  in  the  channel  and  is  given  by 

d(M Yk)  d(JkWk )  .  ...  . 

V  dt  L  =  -  Qy  U  +  skWkAs.  (8) 

Here,  0  is  the  porosity,  sk  is  the  heterogeneous  molar  production 
rate  of  the  chemical  species  k,  y  is  the  independent  spatial  variable 
along  the  thickness,  and  As  is  the  specific  catalyst  area  available  for 
surface  reactions.  Total  density  of  the  fluid  within  the  porous 
structure  can  be  calculated  from 


The  first  term  on  the  right-hand  side  of  Eq.  (10)  represents  the 
diffusive  flux  and  the  second  term  represents  the  viscous  flux. 
Dj}GM  is  defined  as  DGM  diffusion  coefficients  given  as  [46] 

ddgm  =  h-\  (11) 

where  the  elements  of  the  H  matrix  are 


hki  = 


Dl  Kn 


j  *k  Ukj 


&kl  +  ( -  1)  rjF- 
ukl 


(12) 


The  permeability  Bg  in  Eq.  (10)  is  given  by  Kozeny— Carman 
relationship  [50], 


72t(1  -  cp )2  ' 


(13) 


Here,  dp  is  the  particle  diameter  and  t  is  the  tortuosity. 
The  Knudsen  diffusion  coefficient  DfI<n.Eq.  (10)  is  given  by 


% n 


rdp  /& RT 
t  3  y  ttWV 


(14) 


The  global  charge  transfer  reaction  at  the  three-phase  interface 
can  be  written  as: 


H2  +  02-«H20  +  2e-.  (15) 

The  anode  chemistry  model  adopted  here  emphasizes  the 
potential  of  elementary  heterogeneous  mechanism  for  the  steam 
reforming  of  methane  on  Ni  catalysts.  The  modeling  framework  of 
heterogeneous  chemistry  is  described  in  our  previous  reports 
[46,49,51  ].  The  computational  procedure  used  to  solve  the  system 
of  coupled  non-linear  equations  follows  a  space  marching  algo¬ 
rithm,  so  that,  at  each  axial  position,  the  transient  systems  of 
equations  are  solved  until  a  steady  state  solution  is  obtained.  The 
initial  condition  at  each  axial  position  assumes  the  converged 
solution  from  the  previous  finite  volume  cell.  More  discussions 
about  the  procedure  are  given  elsewhere  [52].  Eqs.  (1)— (4)  and  (8) 
form  a  set  of  differential  algebraic  equations  at  each  time  step, 
which  are  solved  after  discretization,  using  the  solver  LIMEX  [55]. 


2.2.  Stacks 


The  timescales  of  various  processes  such  as  kinetics,  diffusion, 
and  heat  transfer  occurring  in  an  SOFC  stack  are  different  from  each 
other,  and  the  heat  transfer  process  has  a  higher  time  constant 
compared  to  the  rest  of  the  processes.  Therefore,  it  is  quite  rational 
to  assume  all  processes  but  heat  transfer  to  be  in  steady  state.  In  the 
method  adopted  here,  the  solid  phase  temperature  is  decoupled 
from  the  fluid  phase  to  develop  the  transient  stack  model,  which 
solves  the  transient  two-  or  three-dimensional  heat  conduction 
equation.  While  solving  the  heat  conduction  equation,  the  stack  is 
assumed  to  be  a  porous  media  consisting  of  straight  channels. 


/<= i  /<= i 


dt 


(9) 


In  the  above  equations,  the  fluxes  Jk  are  evaluated  using  the 
DGM.  According  to  DGM,  the  net  species  molar  flux  is  given  by  [46] 


r  k. 


jk  =  - 


ED«GMv[x']+  E 

L/=i  \/=i 


Dh<n  J 


Vp 


(10) 


_a_ 

dXj 


+  <?> 


(16) 


where  t  is  the  time,  T  is  the  temperature,  p  is  the  density,  Cp  is  the 
heat  capacity,  Xy  is  the  tensor  of  heat  conductivity,  q  is  the  heat 
source  term  from  the  interaction  with  the  channels.  Time  integra¬ 
tion  is  carried  out  for  the  solid  phase  heat  balance.  At  every  axial 
position,  solid  phase  temperature  is  obtained  and  subsequently 
used  to  calculate  the  reaction  rates  and  gas  phase  temperature.  The 
MEA  is  collectively  considered  as  porous  media  by  the  solid 
structure. 
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2.3.  Boundary  conditions 

Solving  the  partial  differential  equations  requires  the  specifi¬ 
cation  of  boundary  conditions  at  the  borders  of  the  simulation 
domain.  Depending  on  the  operating  conditions,  Dirichlet 
boundary  conditions 


T|wan  =  const,  (17) 

or  Neumann  conditions 


k|l 

On 


=  W, 


wall 


(18) 


may  be  used.  In  the  above  equations,  the  heat  flux  W  in  the  normal 
direction  of  the  border  n  may  consist  of  a  linear  heat  condition 
term,  a  bi-quadratic  radiation  term  and  a  constant  term 


V  =  h(T-  Tsurr)  +  EO (t4  -  7^urr)  +  W const •  (19) 

Here,  a  is  the  Stefan-Boltzmann  constant.  The  heat  transfer 
coefficient  h,  the  emissivity  r,  and  the  temperature  of  the 
surroundings  Tsurr  have  to  be  specified  as  model  parameters. 
Nevertheless,  the  stack  model  may  easily  be  coupled  with  an 
arbitrary  combination  of  user-defined  Dirichlet  or  Neumann 
boundary  conditions.  The  initial  condition  used  in  this  case  refers  to 
the  initial  temperature  of  the  solid  structure/stack.  Physically,  it  is 
the  operating  temperature  to  which  the  stack  is  heated  up  before 
use,  since  SOFCs  operate  at  a  range  of  temperatures  from  600  to 
1000  °C. 


2.4.  Heat  source  term 

The  heat  source  term  is  derived  from  the  simulation  of  indi¬ 
vidual  cells.  For  systems  with  nearly  constant  pressure,  the  energy 
conservation  can  be  expressed  in  terms  of  enthalpy.  Constant 
pressure  is  assumed  while  modeling  the  individual  channels.  Then, 
the  enthalpy  flux  in  the  channel  HChannei  can  only  be  changed  along 
the  channel  axis  by  heat  exchange  with  the  stack.  If  the  channel 
density  is  y  (channels  per  unit  area  of  the  cross-section),  the  source 
term  can  be  expressed  as 

(?  =  _Y.9H^nnel  +  Qohmic  (20) 

Here,  Qohmic  is  the  heat  release  due  to  ohmic  heating  within  the 
electrolyte.  By  inclusion  of  enthalpy  of  formation  into  its  definition, 
^channel  does  not  change  due  to  electrochemical  oxidation. 

3.  Solution  methods 

3.1.  Discretization  of  the  PDE 

The  stack  is  assumed  to  be  a  continuum,  whose  transient 
temperature  field  is  described  by  the  heat  conduction  equation  Eq. 
(16).  This  differential  equation  shall  be  solved  for  the  two-  and  three- 
dimensional  cases.  We  assume  that  all  channels  of  the  unit  cells 
within  the  stack  are  parallel,  denoted  by  the  direction  of  the  x-axis. 

The  two-  or  three-dimensional  fields  need  to  be  discretized  in 
order  to  solve  the  differential  equations.  Here,  we  chose  a  finite 
volume  method.  For  the  two-dimensional  case,  we  always  use  an 
orthogonal  grid.  In  case  of  three  dimensions,  the  grid  in  the  yz- 
plane  can  be  unstructured  and  consisting  of  arbitrarily  shaped 
triangles  and  quadrilaterals.  This  grid  is  extended  orthogonally  into 
the  x-direction.  Thus,  the  finite  volume  elements  are  three-  or  four¬ 
sided  prisms. 


Due  to  the  orthogonality  of  the  grid  with  respect  to  the  direction 
of  the  unit  cell  channels,  it  is  easy  to  distinguish  between  the  axial 
and  radial  direction  for  the  heat  conduction.  In  order  to  account  for 
structure  dependent  differences,  we  applied  separate  heat 
conductivity  coefficients  to  the  axial  and  the  radial  directions,  Aaxiai 
and  Aradiai,  respectively. 

The  heat  source  term  q  is  obtained  from  the  simulations  of  indi¬ 
vidual  cells.  Each  row  of  grid  cells  in  the  axial  direction  defines 
a  temperature  profile  that  is  used  as  the  boundary  condition  for  the 
channels  crossing  these  grid  cells.  Thus,  the  number  of  cells  of  the  grid 
in  the  y-direction  or  in  the  xz-plane  defines  a  maximum  number  of 
different  temperature  profiles  for  the  unit  cell  simulations. 

In  order  to  resolve  the  heat  field,  Eq.  (16)  is  discretized  by  finite 
volume  method.  For  a  partial  differential  equation  (PDE)  of  the 
form 


(21) 


we  integrate  over  a  channel  and  apply  Stokes’  integral  theorem  to 
get 


V  pCpfr  =  £  ^ce  (22) 

faces 

Here,  4>face  is  the  heat  flow  through  a  face  of  the  channel.  The 
channel  itself  shall  be  denoted  by  index  0  and  its  neighbor  by  index 
n.  In  the  axial  direction,  we  can  directly  use  a  finite  difference 
approximation  of  the  temperature  gradient  because  of  the 
orthogonality  of  the  grid 


(23) 

In  radial  direction,  the  gradient  can  be  approximated  by  the 
temperatures  T0,  Tn,  TA,  TB,  as  shown  in  Fig.  1.  The  heat  flow  then 
becomes 


^rad  =  (24) 

with 

(25> 

Vaz/ 

This  expression  can  be  simplified  into 

$AB  =  -k(«(T0  -  Tn)  +  0(Tb  -  Ta)).  (26) 

with 


Fig.  1.  Depiction  of  gradient  approximation  in  radial  direction  by  specified 
temperatures. 
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a 


- >2 

AB 


AB  x  PnP0 1 


(27) 


0  =  AB  P^  •  (28) 

|  AB  x  PnP0 1 

The  parameters  a  and  (3  only  have  to  be  calculated  once  for  the 
grid.  As  a  result  of  finite  volume  approximation,  these  grid-based 
parameters  depend  on  Po,  Pn,  which  refer  to  centers  of  adjacent 
cells,  and  A,  B,  which  are  nodes  of  the  same  cell.  With  them,  the 
right-hand  side  of  the  PDE  is  transformed  into  a  system  of  algebraic 
equations. 


3.2.  Choosing  representative  channels 


*k=(rUi  .-.rUn)-  (29) 

of  the  kth  discretized  channel  to  an  output  vector  containing  the 
heat  source  terms  q  of  the  respective  grid  points.  This  mapping  can 
be  expected  to  be  continuous;  channels  with  similar  input  vectors 
shall  be  similar  in  source  terms.  For  the  identification  of  similar 
channels,  an  agglomerative  cluster  algorithm  is  applied. 

A  weight  wk  can  be  assigned  to  each  vector  xk  that  accounts  for 
the  absolute  number  of  channels  represented  by  xk,  which  is 
proportional  to  the  size  of  the  corresponding  stack  area.  In  addition, 
a  distance  function  d(x\  xJ)  is  necessary,  for  which  a  normalized 
maximum  norm  is  chosen: 

d(2,xJ j  =  max(dr(T^alu,r^1),...,dr(T^alln,r^n)). 

(30) 

with  distance  functions 


Calculating  individual  channel  flow  fields  is  the  most  time 
consuming  step  in  a  stack  simulation.  In  order  to  reduce  the 
computational  cost,  only  representative  channels  are  chosen  for 
detailed  simulation.  Due  to  the  representation  of  the  stack  by 
a  discrete  grid,  there  is  a  unique  channel  boundary  temperature 
profile  for  each  axial  row  of  grid  points.  A  finite  number  of  variables 
(for  instance,  only  temperature  T  is  considered  in  this  case) 
completely  define  the  parameters  for  a  channel  simulation,  i.e.,  the 
calculation  of  source  terms  can  be  viewed  as  a  mapping  from  the 
input  vector 


/  .  A  \v  -  V\ 

dT(rjJ )  =  1  gT  1  (31) 

i.e.,  temperatures  are  normalized  with  respect  to  a  temperature 
difference  5T.  The  agglomerative  cluster  algorithm,  illustrated  in 
Fig.  2,  can  be  sketched  as  follows: 

1.  Find  the  minimum  d *  =  d(x*,x*)  =  minI>J(d(xI,xJ)). 

2.  Stop  if  d*  >  1, 

3.  Let  x*  =  (w*x*  +  w*x*)/(w*  +  w*)  and  w*  =  w*  +  w*. 


Fig.  2.  Agglomerative  cluster  algorithm,  (a)  Given  set  of  data  points,  (b)  Join  two  data  points  into  one  cluster,  (c)  After  the  third  agglomeration  step,  two  clusters  have  been  formed, 
(d)  Final  representation  of  the  data  by  four  clusters. 
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4.  Eliminate  (x*,w*),(x*,w*)  and  replace  by  (x*,w*). 

5.  Go  to  1. 

The  remaining  xk  are  the  input  vectors  for  representative 
channels.  All  input  vectors  in  one  cluster  are  associated  with  the 
same  output  source  terms.  Since  the  channel  equations  are  not 
explicitly  time-dependent,  the  number  of  channel  calculations  can 
be  further  reduced  by  including  the  xk  of  previously  calculated 
channels  into  the  clustering  and  reusing  their  results.  A  schematic 
representation  of  the  coupling  between  the  solid  section  and  the 
individual  channel  is  shown  in  Fig.  3.  The  stack  code  is  available  as 
a  part  of  DETCHEM  software  [48]. 


4.  Results  and  discussion 

For  demonstration  purposes,  we  consider  a  stack  having  an 
active  area  of  10  cm  x  10  cm  and  consisting  of  30  cells  in  co-flow 
configuration.  The  fuel  and  air  channels  are  assumed  to  be  of 
1  mm2  cross-sectional  area.  All  the  physical  parameters  and 
thermal  properties  used  to  describe  the  stack  are  given  in  Table  1. 
The  electrochemical  model  parameters  used  for  the  calculations  are 
shown  in  Table  2.  The  exchange  current  density  formulations  are 
taken  from  Ref.  [49].  A  detailed  discussion  about  the  electro¬ 
chemical  model  used  in  this  study  is  enunciated  in  Ref.  [46,49]. 
However,  the  exchange  current  density  parameters  are  adjusted  to 
produce  ~  0.72  V  at  300  mA  cm-2.  Adiabatic  and  Neumann 
boundary  conditions  are  implemented  for  the  present  calculation. 
It  should  be  noted  that  the  time  constants  pertaining  to  electro¬ 
chemistry,  diffusion,  internal  impedance,  mass  transfer  and 
temperature  dynamics  are  specific  to  the  model  parameters  used  in 
this  study. 

As  described  in  the  previous  section,  either  the  voltage  or  the 
current  can  be  fixed  to  estimate  the  rest  of  the  electrochemical 
parameters.  Therefore,  two  different  transient  responses  following 
load  changes  are  studied  here:  (i)  load  change  at  constant  current 
and  (ii)  load  change  at  constant  voltage.  For  both  cases,  the  stack  is 


Table  1 

Stack  component  parameters  and  properties. 


Physical  Flow  channels 

Length  (m) 

0.1 

Height  (m) 

1  x  10-3 

Width  (m) 

Anode 

1  x  10-3 

Thickness  (pm) 

750 

Porosity  (%) 

35 

Tortuosity 

3.5 

Particle  diameter  (pm) 

2.5 

Pore  diameter  (pm) 

1.0 

Specific  area  (m-1) 

Cathode 

1.080  x  105 

Thickness  (pm) 

30 

Porosity  (%) 

35 

Tortuosity 

3.5 

Particle  diameter  (pm) 

2.5 

Pore  diameter  (pm) 

1.0 

Specific  area  (m-1) 

Electrolyte 

1.080  x  105 

Thickness  (pm) 

15 

Thermal  Solid/stack 

Specific  heat  capacity  (J  kg-1 1<) 

440 

Density  (kgm-3) 

5940 

Thermal  conductivity  (J  m-1  s  K) 
Insulator 

1.86 

Specific  heat  capacity  (J  kg-1 1<) 

1047 

Density  (kgirT3) 

480 

Thermal  conductivity  (J  m-1  s  K) 

0.059 

Thermal  emissivity 

0.09 

Convective  coefficient  insulator/ 
surrounding  (W  m-2 I<-1) 

2 

initially  assumed  to  be  at  800  °C.  Inlet  fuel  and  air  are  assumed  to 

be  at  800  °C  and  750  °C,  respectively. 

4.1.  Constant  current 

4.1.1.  Adiabatic  condition 

Fig.  4  displays  the  transient  response  in  the  stack  temperature 
for  cells  operating  at  constant  current.  From  the  initial  condition, 
the  temperature  distribution  in  the  stack  reaches  steady  state  after 
~  14.5  min.  At  t  =  80  min,  a  step  change  in  current  from  0.3  A  cm-2 
to  0.4  A  cm-2  is  introduced.  As  expected,  the  stack  temperature 
increases  with  the  step  increase  in  current  and  reaches  steady  state 
after  another  ~  17  min.  These  response  times  correspond  well  with 
the  values  reported  by  Ivers-Tiffee  et  al.  [33].  The  corresponding 
changes  in  cell  voltage  and  reversible  voltage  are  shown  in  Fig.  5. 
The  reversible  potential  decreases  following  the  step  change  in 
current  due  to  the  increase  in  temperature.  However,  subsequent  to 

the  step  change  in  current  from  300  mAcm  2 

to  400  mA/cm2,  the 

cell  voltage  first  decreases  before  reaching  steady  state.  The  initial 
decrease  in  cell  potential  occurs  due  to  the  following  reasons. 

•  An  increase  in  current  leads  to  an  increase  in  overpotential 

losses,  because  the  stack  temperature  at  the  time  of  step 

Table  2 

Electrochemical  model/input  parameters. 

Anode  asymmetry  factor  (/3a) 

0.5 

Cathode  asymmetry  factor  (/?c) 

0.5 

Exchange  current  density  parameters 

Pre-exponential  for  H2  oxidation  (/<H2)  (A  cm-2) 

1.07  X  104 

Pre-exponential  for  02  reduction  ( k02 )  (A  cm-2) 

5.19  x  103 

Activation  energy  for  H2  oxidation  (£H2)  (J  mol-1) 

87.8  x  103 

Activation  energy  for  02  reduction  (E02)  (J  mol-1) 

88.6  x  103 

Fig.  3.  Illustration  of  coupling  between  the  solid  structure  and  the  channel. 
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Fig.  4.  Maximum,  minimum,  and  average  temperature  in  the  stack  for  step  change  in 
current  from  0.3  A  cm  2  to  0.4  A  cm-2. 


change  still  corresponds  to  the  steady  state  temperature  at 
300  mA  cm-2. 

•  The  activation  losses  gradually  decrease  as  the  temperature 
approaches  the  new  steady  state. 

In  short,  the  under-shoot  in  cell  potential  originates  from  the 
over-shoot  in  activation  losses  following  the  step  change.  The  cor¬ 
responding  transient  response  of  activation  losses  at  the  anode  and 
cathode  are  displayed  in  Figs.  6  and  7,  respectively. 

In  Fig.  8,  we  see  the  response  of  the  stack  temperature  to  step 
changes  in  current  density  from  0.3  A  cm-2  and  0.5  A  cm-2  to  a  new 
steady  state  value  of  0.4  A  crrT2.  The  response  times  of  the  system, 
after  the  respective  step  changes  are  introduced,  are  the  same  and 
equate  to  ~  17  min.  Thus,  we  observe  that  the  response  times 
remain  independent  of  the  magnitude  of  load  introduced.  This 
phenomenon  was  also  observed  by  Achenbach  [37]. 
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Fig.  6.  Anode  activation  overpotential  following  a  step  change  in  current  from 
0.3  A  cm-2  to  0.4  A  cm-2. 
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Fig.  7.  Cathode  activation  overpotential  following  a  step  change  in  current  from 
0.3  A  cm-2  to  0.4  A  cm-2. 


Fig.  5.  Cell  potential  and  reversible  potential  following  a  step  change  in  current  from 
0.3  A  cm-2  to  0.4  A  cm-2. 


Fig.  8.  Transient  temperature  distribution  during  load  change  from  j  =  0.3  or  j  =  0.5  to 
j  =  0.4  A  cm-2. 
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Fig.  9.  Temperature  distribution  in  the  3D  stack  (solid  phase)  at  different  surface 
planes,  along  and  across  the  flow  direction,  at  t  =  60  min  during  adiabatic  operation. 


The  temperature  distribution  in  the  stack  at  different  planes 
along  and  across  the  flow  direction,  after  the  first  steady  state,  is 
displayed  in  Fig.  9.  As  expected,  the  temperatures  in  the  stack 
increase  in  a  monotonic  manner  along  the  flow  direction  due  to  the 
combined  effect  of  exothermic  cell  reactions  and  the  current 
consumed  by  the  internal  cell  resistance.  However,  uniform  profiles 
are  observed  across  the  direction  of  flow  because  of  the  relatively 
thin  PEN  structure  and  the  radial  thermal  conductivity  of  the  solid. 
Also,  the  insulator  seems  to  have  minimal  effect  on  the  internal 
temperature  distribution  of  the  system.  It  should  be  noticed  that 
the  model  has  a  limitation  that  does  not  account  for  individual 
mechanical  components  in  the  stack.  It  is  a  lumped  parameter 
model  as  far  as  individual  mechanical  components  are  concerned, 
and  therefore,  it  is  assumed  that  the  entire  stack  has  the  same 
thermal  properties. 


Table  3 

Solid  structure/stack  material  properties. 

Material  Thermal  conductivity  (A  —  W  mK  1 

)  Thermal  diffusivity  ( a  -mm2  s  a) 

(a)  1.86 

0.716 

(b)  5.86 

2.466 

(c)  27 

6.923 

0.035 

0.033 

0.031 

%  0.029 

G 

V  0.027 

o 

I  0.025 


Q> 
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Fig.  11.  Ohmic  polarization  following  a  step  change  in  current  from  0.3  A  cm  2  to 
0.4  A  cm-2. 


During  the  calculations,  we  observed  that  the  time  constants  for 
the  dynamic  responses  are  strongly  dependent  on  thermal  diffu- 
sivity  (A/pCp).  Lower  thermal  diffusivity  led  to  a  faster  response  or  to 
a  shorter  time  constant,  which  can  be  attributed  to  mass  transfer 


0.08 


0.06 

Jj 

N 

0.04 


XM 


Temperature  (K) 

11200 
1150 

m  1100 
1050 
1000 
950 
900 
850 

U800 

750 

700 


Fig.  12.  Temperature  distribution  in  the  3D  stack  (solid  phase)  at  different  surface 
planes,  along  and  across  the  flow  direction,  at  t  =  60  min  during  computation  with  an 
imposed  Neumann  condition. 


V.  Menon  et  al.  /  Journal  of  Power  Sources  214  (2012)  227-238 


235 


Fig.  13.  (a)  Maximum,  minimum,  and  average  temperature  in  the  stack  for  step  change  in  voltage  from  0.7  V  to  0.8  V.  (b)  Average  current  density  following  a  step  change  in  voltage 
from  0.7  V  to  0.8  V.  (c)  Reversible  cell  potential  following  a  step  change  in  voltage  from  0.7  V  to  0.8  V.  (d)  Anode  activation  overpotential  following  a  step  change  in  voltage  from 
0.7  V  to  0.8  V.  (e)  Cathode  activation  overpotential  following  a  step  change  in  voltage  from  0.7  V  to  0.8  V. 


Fig.  14.  Internal  steady  state  temperature  distribution  of  the  stack  before  step  change  in  voltage  from  0.7  V  to  0.8  V. 
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Fig.  15.  Internal  steady  state  temperature  distribution  of  the  stack  after  step  change  in  voltage  from  0.7  V  to  0.8  V. 


dynamics.  This  was  also  observed  by  Achenbach  [37].  It  might  be 
possible  that  Li  et  al.  [32]  used  high  thermal  diffusivity  for  their 
model  calculations,  which  led  to  very  high  response  times.  This 
phenomenon  is  depicted  in  Fig.  10,  with  the  corresponding  solid 
stack  structure  properties  listed  in  Table  3.  A  number  of  other  stack 
operating  conditions,  such  as  convection  and  radiation  losses,  and 
the  effect  of  re-circulation  zones  at  the  inlet  and  outlet  faces/ 
manifolds  are  not  discussed  here. 

4.1.2.  Neumann  condition 

The  Neumann  condition,  as  described  in  Table  1,  is  applied  at  the 
four  faces  of  the  stack,  barring  the  inlet  and  outlet  face  which  are 
kept  adiabatic.  As  shown  in  Fig.  4,  the  temperature  distribution 
reaches  steady  state,  from  the  initial  value,  at  -15  min.  A  step 
change  in  current  from  0.3  A  cm-2  to  0.4  A  cm-2  at  t  =  80min  is 
introduced,  due  to  which  the  temperature  distribution  in  the  stack 
reaches  a  new  steady  state  value  in  -  18.5  min.  This  is  due  to  the 
dependence  of  response  time  on  operating  conditions. 

The  heat  loss  at  the  surface  of  the  stack  results  in  a  decrease  in 
the  magnitude  of  its  temperature  distribution,  which  leads  to  the 
production  of  relatively  low  amount  of  waste  heat.  This  causes  the 
average  cell  temperature  to  decrease,  which  manifests  itself  as  an 
increase  in  ohmic  polarization,  resulting  in  a  lower  operating  cell 
voltage  as  described  in  Fig.  11.  The  reversible  potential,  which  is 
obtained  from  the  Nernst  equation,  is  higher  due  to  the  lower 
amount  of  heat  generated  by  the  cell  (TA5)  that  operates  reversibly. 
This  increases  the  amount  of  free  energy  available  to  the  system. 
This  phenomenon  is  shown  in  Fig.  5. 


The  internal  temperature  distribution  in  the  stack,  at  different 
planes  after  the  attainment  of  the  first  steady  state,  is  shown  in 
Fig.  12.  This  is  a  pre-dominant  factor  in  determining  any  induced 
thermal  stresses  in  the  materials  positioned  at  different  locations  of 
the  system.  Furthermore,  as  a  result  of  heat  loss/flux  at  the 
boundaries,  the  distinct  temperature  contours  are  formed  due  to 
the  increase  in  temperature  gradients  in  the  Yand  Z  directions.  This 
shows  the  importance  of  3D  modeling  for  determining  how 
temperature  dynamics  play  an  important  role  in  the  transient 
responses  of  a  SOFC  system. 

4.2.  Constant  voltage 

4.2.1.  Adiabatic  condition 

Fig.  13a  describes  the  transient  response  in  temperature  for  cells 
operating  at  constant  voltage.  Flere,  the  initial  steady  state  is  ach¬ 
ieved  in  - 15  min,  and  at  t  =  80  min,  a  step  change  in  voltage  from 
0.7  V  to  0.8  V  is  introduced.  The  increase  in  cell  potential  decreases 
the  average  current  from  the  stack,  and  as  expected,  the  stack 
temperature  reaches  a  new  lower  steady  state  value  in  about 
— 18  min.  The  corresponding  transient  responses  in  current  and 
reversible  potential  are  displayed  in  Fig.  13b  and  c,  respectively.  As 
usual,  the  decrease  in  temperature  following  the  step  increase  in 
voltage  leads  to  an  increase  in  the  reversible  potential.  The  current 
density  drops  from  —340  mA  cm-2  following  the  step  change  in 
voltage,  and  then  gradually  assumes  a  new  steady  state  value  of 
-150  mAcnrr2.  Nevertheless,  the  step  change  in  voltage  leads  to 
an  under-shoot  in  activation  losses.  At  the  time  of  step  change,  the 
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temperature  in  the  stack  corresponds  to  that  of  ~ 340  mA  cm-2; 
however,  the  drop  in  current  density  at  higher  stack  temperature 
leads  to  lower  activation  losses,  which  in  turn  corresponds  to  the 
point  of  under-shoot.  Nonetheless,  as  the  stack  temperature 
assumes  a  new  lower  steady  state  value,  the  activation  losses 
increase  and  asymptotically  reach  a  new  steady  state  value.  The 
anode  and  cathode  activation  losses  are  displayed  in  Fig.  13d  and  e, 
respectively. 

4.2.2.  Neumann  condition 

The  thermal  properties  of  the  insulator  and  the  Neumann 
condition  at  the  stack  surface  are  described  in  [54].  Fig.  13a  shows 
that  the  time  taken  to  reach  steady  state  from  the  initial  condition 
is  ~17  min.  At  t=  80  min,  a  step  change  in  voltage  from  0.7  V  to 
0.8  V  is  introduced,  causing  the  stack  temperature  to  attain  a  new 
steady  state  value  in  about  ~  17  min.  The  ohmic  polarization 
increases  as  the  operating  temperature  decreases,  which  results  in 
the  decrease  of  current  density  at  a  given  cell  voltage,  as  shown  in 
Fig.  13b.  Consequently,  the  reversible  potential  increases  as  shown 
in  Fig.  13c. 

The  3D  contour  plots  of  the  cell  stack  describe  the  steady  state 
internal  temperature  distributions  before  and  after  the  step  change 
in  voltage  is  introduced,  which  are  shown  in  Figs.  14  and  15, 
respectively.  These  figures  represent  the  temperature  profiles  of  the 
cell  located  at  the  center  of  the  stack.  A  comparison  with  the 
adiabatic  case  provides  a  comprehensible  idea  about  the  differ¬ 
ences  in  the  magnitudes  of  temperature  gradients  for  the  same 
model  input  parameters.  It  can  be  seen  that  the  adiabatic  condition 
results  in  a  smaller  temperature  gradient  at  the  boundaries 
compared  to  the  case  involving  heat  loss  at  the  boundaries.  The 
resulting  temperature  distribution  arises  due  to  the  net  effect  of 
heat  absorbed  or  released  as  a  result  of  heterogeneous  chemical 
reactions,  heat  release  due  to  the  electrochemical  reactions  at  the 
TPB,  resistive  heating  within  the  electrolyte,  convective  heat 
transfer  to  and  from  the  channels,  radiation  heat  transfer  with  the 
interconnects,  and  also  heat  loss  at  the  stack  boundaries  in  the  case 
of  an  imposed  Neumann  condition.  The  temperature  decline  at  the 
inlet  is  attributed  to  endothermic  reforming  reactions.  It  is  evident 
that  the  temperature  increases  in  the  axial  direction,  due  to 
exothermic  cell  reactions,  along  with  the  air  and  fuel  channel 
temperature.  This  in  turn  reflects  in  the  gradual  increase  of  insu¬ 
lator  temperature  from  the  inlet  to  the  outlet  face.  Each  cell’s 
temperature  distribution  also  depends  on  its  position  within  the 
stack.  This  affects  the  dynamics  of  almost  all  other  mechanisms. 
During  our  simulations,  when  considering  different  materials  for 
insulation  at  the  stack  surface,  we  observed  that  the  choice  of 
insulation  layer  was  mainly  responsible  for  controlling  the  surface 
temperature  of  the  system  that  determined  the  amount  of  heat 
available  for  transfer  to  the  surroundings. 

5.  Conclusion 

We  have  developed  a  novel  numerical  algorithm  for  the 
modeling  of  3D/2D  SOFC  stacks.  The  model  is  demonstrated  with 
a  reformate  fuel  composition  for  a  co-flow  cell  configuration. 
Instead  of  simulating  every  individual  cell  in  the  stack,  only 
representative  ones  based  on  the  cluster  algorithm  are  chosen  for 
detailed  simulation.  This  approach  dramatically  reduces  the 
computation  time  required  for  stack  simulation.  Furthermore,  we 
also  studied  the  dynamics  of  the  stack  during  two  different  tran¬ 
sient  responses  to  load  changes:  constant  current  and  constant 
voltage.  It  was  observed  that  the  adiabatic/Neumann  cases  at  the 
borders  of  the  simulation  domain  yielded  almost  the  same 
response  times  between  the  two  steady  states,  for  both  types  of 
load  changes.  The  limitation  of  the  model  is  the  higher  level  of 


abstraction  related  to  the  individual  mechanical  components  of  the 
stack.  Nevertheless,  the  novel  model  and  computational  tool  can 
assist  in  a  better  understanding  of  heat  balances  in  SOFC  stacks, 
even  when  operated  with  other  fuels  than  pure  hydrogen. 
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Nomenclature 

Ac\  area  of  cross-section  of  flow  channels  (m2) 

As:  specific  area  (m-1) 

Bg\  permeability  (m2) 

Cp:  specific  heat  (J  kg-1  I<-1) 
dp:  particle  diameter  (m) 

D:  diffusivity  (m2  s-1) 

Del:  effective  binary  diffusion  (m2  s-1) 

Dh:  hydraulic  diameter  (m) 

Dft Kn:  effective  Knudsen  diffusion  (m2  s-1) 


ECeii-  cell  voltage  (V) 

F:  Faraday  number  (C  mol-1) 

Gz:  Graetz  number 

h:  heat  transfer  coefficient  (J  m-2  K-1  s-1);  specific  enthalpy  (J  kg-1) 
Hc:  channel  height  (m) 

AH:  enthalpy  of  formation  (J  mol-1) 

H:  enthalpy  flux  (J  s-1) 
i:  current  density  (A  cm-2) 

Ji<:  species  flux  (mol  m-2  s-1) 

I<g :  number  of  gas-phase  species 
n:  number  of  charge  transferred 
Nu:  Nusselt  number 
p:  pressure  (Pa) 

Pe:  MEA  perimeter  (m) 

Pr:  Prandtl  number 

Qohmic-  ohmic  heat  source  (J  m-3  s-1) 

R:  gas  constant  (J  mol-1  K-1) 

Re:  Reynolds  number 

s:  molar  production  rate  (mol  m-2  s-1,  mol  m-3  s-1) 

AS:  entropy  change  (J  K-1) 
t:  time  (s) 

T:  temperature  (K) 
u:  velocity  (ms-1) 

V:  volume  (m3) 

V14:  molecular  weight  of  fcth  species  (kg  mol-1) 

[X]:  concentration  (mol  m-3) 
x,  y,  z:  coordinate  direction  (m) 

X:  mole  fraction 
Y:  mass  fraction 


Greek  symbols 
y:  channel  density  (m-2) 

5:  Kronecker  delta  symbol 
r:  emissivity 

X:  thermal  conductivity  tensor  (J  m-1  s-1  K-1) 

fi:  viscosity  (kg  m-1  s-1) 

p:  density  (kgm-3) 

ae:  conductivity  (Scm-1) 

a:  Stefan— Boltzmann  constant  (W  m-2  K-4) 

t:  tortuosity 

0:  porosity 

W:  heat  flux  (J  m-2  K-1) 

Subscripts 
c:  channel 

e:  electrolyte/electrode 
f:  fluid 

I:  interconnect 
k:  species  index 


